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Abstract 

Lagrangian coherent structures (LCS) in fluid flows appear as co-dimension 
one ridges of the finite time Lyapunov exponent (FTLE) field. In three- 
dimensions this means two-dimensional ridges. A fast algorithm is presented 
here to locate and extract such ridge surfaces while avoiding unnecessary 
computations away from the LCS. This algorithm reduces the order of the 
computational complexity from 0{l/dx^) to about by eliminating 

computations over most of the three dimensional domain and computing 
the FTLE only near the two-dimensional ridge surfaces. The algorithm is 
grid based and proofs of error bounds for ridge locations are included. The 
algorithm performance and error bounds are verified in several examples. 
The algorithm offers significant advantages in computational cost as well as 
later data analysis. 
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1. Introduction 

Lagrangian coherent structures (LCS) have seen increasingly popular use 
for visualizing and quantifying fluid structures and transport and mixing 
behavior in fluid flows. LCS were proposed by Haller and Yuan as the locally 
most repelling or attracting material lines in a flow [1]. Shadden et al. |2] 
later defined LCS as ridges of the finite time Lyapunov exponent (FTLE) 



field. Throughout this paper we will use this definition of LCS (ridges of the 
FTLE field), but other definitions of LCS will be discussed in the conclusions. 
LCS are known to have many useful properties such as denoting barriers 
to transport and repelling or attracting material surfaces [2]. Additionally, 
they form unambiguous boundaries to well known coherent structures such 
as vortices [31 E] and are relatively insensitive to small errors [3] . 

The FTLE field is typically computed numerically by integrating particle 
trajectories through the fiow to approximate the fiow map and then using 
finite differences to approximate the Jacobian of the flow map. The flow map 
which maps particles from their initial position at time to to a flnal position 
at time + T, ^, is given by 



rto+T 

$^»+^(x) = x(to) + / v(x(t))dt (1) 

J to 



and may be computed from analytical, experimental, or numerical velocity 
data. If the velocity data is not analytically deflned it is usually necessary 
to read the velocity from data flies and perform interpolations in space and 
time. The Jacobian of $ is used to compute the deformation tensor. A, 
which contains information about the stretching in the flow. 



dx y \ dx y 

Finally, the largest eigenvalue of A is used to deflne the flnite time Lyapunov 
exponent, a, 

rj;(x) = -^lnv/A™,,(A). (3) 



cr. 



T 



Large values of the FTLE correspond to large amounts of stretching in the 
flow. Additionally, one may compute the flow map either forward or back- 
ward in time (positive or negative T). Ridges in the FTLE fleld are then ex- 
pected to correspond to either the locally most attracting or repelling lines in 
the flow. This is not strictly true since it is also possible for high shear regions 
to exhibit large FTLE values, but experience has shown that FTLE ridges are 
often sufficient to learn about the underlying flow structure. More advanced 
techniques may be used to ensure that shear structures are not selected or 
that the resulting LCS are exact barriers to transport (e.g. Haller, [6]). 

One of the largest hurdles to more widespread use of LCS techniques is the 
large time required to compute the FTLE. For typical fluid flows, computing 
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the FTLE field requires advecting large numbers of particles through the 
flow. This must be done for each time step at which the FTLE held is 
desired. These particle advections dominate the total computational cost of 
any LCS algorithm. 

More efficient algorithms are needed to address the large cost associated 
with LCS computations. Several past papers have attempted to address this 
problem with adaptive mesh refinement (AMR) algorithms that refine the 
computational mesh near the LCS [ZIIH]- This is effective, but still results in 
computing many FTLE values away from the LCS. 

Additional attempts have been made to reuse computations to compute 
the LCS at subsequent time steps [H]. Since the LCS are often desired at 
many different times and the time step may be smaller than the integration 
time used this leads to overlapping integrations. For example, if the LCS 
are desired at times t = {0.0, 0.1, 0.2, ...10.0} and an integration time of 
T = 1.0 is to be used, the standard approach would be to compute a series 
of flow maps {$o!o°' "^oV' ^o^i^y ■ ■ ■} each time the FTLE field is required. 
However, it is instead possible to compute the fiow maps {^q I, $o;i, ^o'^, ■ ■ ■} 
and use the fact that $J°o° = "^gV ° ° ^l'^ ' ' ' "^ao- This composition of 
fiow maps technique has the potential for large efficiency gains, but only if 
there is a significant overlap in the integration times and many time steps 
are desired. The method also comes at the cost of greatly increased memory 
usage to store all the necessary fiow maps. 

Another recent paper has reformulated the FTLE problem as an Eule- 
rian level set problem involving the solution of a Liouville equation [TO] . 
This allows the use of any previously developed high order accurate schemes 
for the resulting hyperbolic PDE's. Although there are many existing tech- 
niques for solving such hyperbolic systems, the time step used in the Eulerian 
method must obey a CFL condition to ensure stability while the time step in 
more commonly used Lagrangian techniques is typically determined by the 
required accuracy. This typically means that the Lagrangian techniques are 
faster. 

We propose another alternative to any of these approaches: detect and 
track the ridges on the fiy. A recent publication discusses a gridless ridge 
tracking algorithm for computing the FTLE ridges in 2D fiows [11]. How- 
ever, despite it's effectiveness in 2D (providing speedups of up to 80x), that 
technique is not accompanied by proofs of convergence or error bounds on 
the results and does not readily generalize to higher dimensions. This is 
because one- dimensional ridges in a 2D flow may be represented as a sim- 
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pie curve and need only be tracked in two directions, but in three or more 
dimensions, the ridges are at least two-dimensional surfaces which require a 
more sophisticated representation. 

In this paper, we present a fast algorithm for efficiently computing LCS 
in three-dimensional flows. In fact, most of the algorithm may be used for 
n > dimensions, but the surface triangulation used here is specific to 2D 
surfaces in a 3D space. All computations are performed on a predetermined 
orthogonal grid to simplify implementation and surface triangulation. A 
grid-less algorithm requires significant amounts of time to generate surface 
meshes, but by using a fixed grid, we are able to efficiently generate a surface 
triangulation via a lookup table similar to the marching cubes algorithm that 
is used for computing isosurfaces [12] • FTLE ridges are initially detected by 
computing the FTLE values on a series of lines across the domain. Local 
maxima along these lines occur at the FTLE ridge crossings. Once a series 
of points on the ridges have been detected, nearby points are tested to see if 
they are also on the ridge. This process is repeated to track the ridges through 
the entire domain. By performing computations only near the LCS surfaces, 
the order of the algorithm is reduced from 0{l/dx^) to about 

We present the results from several test cases, including a time dependent 
double gyre, Arnold-Beltr ami- Childress (ABC) fiow, and a swimming jelly- 
fish. These results establish the computational order of the algorithm. We 
also investigate the computation cost for an FTLE field where the surface 
area of the LCS is known a priori and find that for a fixed grid spacing, the 
computational time is C + 0{Alcs) where C is a constant initialization cost 
and Alcs is the area of the LCS surfaces. Finally, the ridge tracking algo- 
rithm offers several other advantages beyond the savings in computational 
time. Since the LCS surfaces are directly computed, visualization of the LCS 
is simplified 

2. A surface tracking algorithm 

The computational savings seen by using a ridge tracking algorithm come 
from avoiding unnecessary computations away from the FTLE ridges. We 
separate this process into three steps: detecting initial points on the ridge 
surfaces, tracking the ridges through space, and triangulating the ridge points 
into a ridge surface. The initial ridge detection is handled by detecting where 
lines through the domain cross the ridge surfaces. The ridges are then iter- 
atively grown by searching for nearby ridge points until no new ridge points 
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Figure 1: Initial ridge detection is handled by looking for local FTLE maxima along a few 
lines through the domain. 

are found. Finally, the use of a gridded coordinate system throughout this 
process allows the use of a lookup table to efficiently generate a triangulation 
of the resulting LCS surfaces. 

To begin, we need a ridge definition. Shadden et al. [2] offers two defi- 
nitions of ridges in two dimensions. Here, we extend the concept of second 
derivative ridges to n dimensions: 

Definition Ridge: A ridge of a function F is a co-dimension one surface 
S satisfying 

1. The vectors n ■ V-F = for all points on S where n is a unit vector 
normal to S. 

2. ri^Hn = min||u||=i (u^Huj < for all points on S where H is the 
Hessian matrix associated with F. 

2.1. Initial ridge detection 

If a hiker walks in a straight line, constantly monitoring his altitude, he 
will reach a locally maximum altitude upon crossing a ridge in the terrain. 
Similarly, if the FTLE is known along a line through a three-dimensional 
domain, local maxima along the line occur where the line crosses FTLE 
ridges. We restrict our computations to a fixed grid with spacing dx and 
initially detect the FTLE ridges by computing the FTLE values along a set 
of lines that cross the domain as seen in Fig. [T] The number and spacing 
of the lines is determined by the user and depends on the expected spatial 
extend of the LCS in the fiow. If an LCS surface does not intersect any of 
the lines used and is isolated from other LCS it may be missed entirely. 
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Once the FTLE values along these lines are known, we look for local 
maxima on each line and define grid ridge points as follows: 

Definition Grid Ridge Point: Given an orthogonal grid in with coordi- 
nate directions : z G {1, .■■,n}; a grid point xq is a grid ridge point of a 
function F if -F(xo) > -F(xo ± e^) for at least one i G {1, ...,n}. 

The local maxima on each line are grid ridge points and we begin tracking 
the FTLE ridges from these points. In the limit as grid spacing goes to zero, 
the grid ridge points converge to coordinate local maxima defined as: 

Definition Coordinate local maximum: A point xq is a coordinate local 
maximum of the function F with respect to coordinate direction ej if there 
is a value e > such that -F(xo) > -F(xo + 5 ■ e^) for all 5 < e. 

This property is proven below in Section |3} 

It is also desirable to set a threshold for the FTLE ridge values at this 
time. Only detecting ridges with FTLE values above some threshold ensures 
that only the strongest LCS are revealed. This is often done by restricting 
the ridges to have FTLE values above a percentage of the maximum FTLE 
values that are detected. Typically 50 — 70% of the maximum value is an 
adequate choice. Although the grid ridge points may exhibit false positives 
in the sense that such points may not necessarily correspond to FTLE ridges 
as defined by Shadden et al. [2], experience has shown that this happens 
infrequently and typically does not alter the topology of the resulting ridge 
surfaces. 

2.2. Ridge tracking 

The heart of this algorithm lies in tracking the ridges outward from the 
initially detected grid ridge points. A schematic of this process for a two 
dimensional example is shown in Fig. [2j The neighbors of each initially 
detected grid ridge point are checked to see if any meet the criteria to be 
grid ridge points (Fig. ^p). In 3D, points in the grid are indexed by (i, j, k) 
for the x, and z directions. For each newly detected grid ridge point, 
(z, j, /c), the neighboring points in [i — 1, z + 1] x [j — 1, j + 1] x [fc — 1, + 1] 
are checked to see if any are grid ridge points. If new grid ridge points are 
found that lie above the FTLE threshold the neighbors of those points are 
checked. This process is repeated until no new grid ridge points are found. 

It is important to note that this algorithm detects grid ridge points as 
defined above, and therefore detects coordinate local maxima (also defined 
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(a) Initial grid ridge points 
detection step, green circles are 
grid ridge points. 



(b) The next set of points to check 
for ridges (black o's) and previously 
checked points (gray o's). 





(c) New grid ridges (blue), new 

points to check (black) and 
previously checked points (gray). 



(d) The resulting LCS. 



Figure 2: The ridge tracking process in 2D. The background color represents the FTLE 
field, +'s mark grid points. The 3D process is analogous, but a surface triangulation is 
used instead of line segments. 
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above). If the ridge height is not constant (i.e. VF 7^ on the ridge) the 
grid ridge points will not converge to the actual ridges. Section [3] provides 
an error bound and corresponding proofs addressing this issue, but the error 
is typically smaller than reasonable grid spacings. 

2.3. Surface triangulation 

At the end of the ridge tracking process, a large list of grid ridge points 
is obtained. These points must be connected into a surface triangulation 
for visualization and further analysis. Initial attempts at developing a 3D 
ridge tracking algorithm revealed that grid-less techniques requiring surface 
meshing where both complex and expensive because of the surface mesh- 
ing process. By using a gridded coordinate system, it is possible to very 
quickly generate a surface triangulation from a lookup table similar to the 
process used in the marching cubes algorithm that is popular for isosurface 
generation [T^ . 

After all the grid ridge points have been found, each cubic element in the 
domain is given an 8-bit number corresponding to the configuration of the 
grid ridge points on the elements 8 vertices. The 256 possible configurations 
are listed in a lookup table which directly converts the 8-bit number to a 
triangulation of the grid ridge points in the element. 

Generating the lookup table is accomplished by using reflections and ro- 
tations to reduce the the 17 canonical cases shown in Fig. [3j Elements 
containing two or fewer grid ridge points contain no LCS surface triangles 
while elements containing 3 or more grid ridge points are added to the surface 
triangulation. 

The lookup table makes the surface triangulation step in this algorithm 
extremely efficient, but at the cost of adding some ambiguity in the precise 
triangulation that should be used for certain ridge point configurations. This 
is a well known problem in isosurface construction and is commonly dealt 
with by making assumptions about the field at the subgrid scale and testing 
additional points within the element. Similar tests may be possible for the 
FLTE ridges, but we have chose to err on the side of adding extra triangles 
to these ambiguous cases surfaces to avoid gaps in the LCS surfaces. 

3. Algorithm properties 

In this section we address the important properties of the ridge track- 
ing algorithm by proving that well defined ridges are detected by the ridge 
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tracking algorithm and the error in the location of the ridges is typically very 
small. The theorems contained below rely on the ridge definition presented 
above in Section [2] A few important properties are summarized here for 
convenience: 

• VF is parallel to the ridge. 

• The eigenvectors {vj} of the Hessian H of F form a complete, orthonor- 
mal basis for the space since H is symmetric. 

• The eigenvector vi associated with the minimum eigenvector Ai of H 
is normal to the ridge. 

To prove that the ridge tracking algorithm accurately detects ridges, we 
first show that grid ridges converge to coordinate local maxima and then show 
that well defined ridges always have a nearby coordinate local maximum and 
find a bound on the distance between a ridge and the nearest coordinate local 
maximum. 

Theorem 3.1. Let F be a C"^ continuous function and let L be a grid line 
parallel to e^. Assume x* on L is a local maximum of F in the e^ direction 
and the second derivative of F in the e^ direction is negative. Then there 
exists a value e such that if the grid spacing dx < e =^ ||x* — x^|| < dx for 
some grid point x^ 

Proof. Since x* is a local maximum in the e^ direction, there is a value ei > 
such that for all |ci| < ei, F|x*+ciei < -F|x*. Since F is continuous and 
d'^F/dei^ < 0, 3 £ e (0,£i) such that 



d^F 



< V |c| < £. (4) 



This implies that dF/dei is monotonically decreasing over the interval c G 
i-e,e). 

Next, choose a grid spacing, dx < e = e/2. This ensures that at least 
five grid points lie in the intervalal (x* — eej, x* + eSi) and at least lie two to 
each side of x*. 

If X* is a grid point, then the neighboring grid points have lower values 
of F since dx < e < ei so this grid point is a grid ridge point and there is a 
grid ridge point less than dx from x*. 
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If X* is not a grid point, label the nearest four grid points xi, x, X3, X4, 
with X* lying between x and X3. Since dF/dei is monotonically decreasing 
over the interval and dF/dei\y^* = 0, F|xi < -F|x and Fjxg > -^"1x4- If 



-^Ix > -^1x35 X is a grid ridge point. If F|x < F\y_^, X3 is a grid ridge point. 
If F|x = -F|x35 both X and X3 are grid ridge points. In any of these cases, 
since x* lies between x and X3, there is a grid ridge point less than dx from 

X*. □ 



Thm. 3.1 proves that as grid spacing goes to zero, there are grid ridge 
points that converge to the coordinate local maxima of F. Therefore the 
ridge tracking algorithm detects coordinate local maxima. We now show 
that every well defined ridge has a nearby coordinate local maximum. 

Theorem 3.2. Given a continuous function F that admits a sufficiently 
sharp second derivative ridge, for every point xq on the ridge, there is a 
nearby point, x* , in the ridge normal direction that is a coordinate local 
maximum, 
ridge than 



The nearest coordinate local maximum is no further from the 
2v/^^||DF(xo)|| 



as long as 



and 



d 



Ivi ■ (D^F)!! < 



A? 



2nv/n^||DF(xo) 



n 



n 



1 



n 



|Ai| 



|DF(xo)|| II vi -0^(011 



where D denotes the gradient operator, Ai (respectively Xn) is the minimum 
(respectively maximum) eigenvalue of the Hessian, D^F(xo), vi and v„ are 
the eigenvectors associated with Ai and A„, and n is the dimension of the 
space. 

Although these criteria may seem restrictive, they are typically easily 
satisfied by well defined ridges as we will see in examples below. Ai is typically 
very large in magnitude (sometimes up to 10^^) while ||DF|| and |A„| are 
typically 0{1). If A„ is negative the third condition is trivially satisfied. 

Proof. By assumption F is continuous and admits a second derivative 
ridge through point xq. Denote the eigenvalues and normalized eigenvectors 
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of the Hessian D^F(xo) as Ai < A2 < ■ ■ ■ < A„ and vi,v,...,v„. By 
definition, vi is normal to the ridge and DF(xo) is parallel to the ridge. We 
also know that {vi, . . . , v„} forms an orthonormal basis for the space since 
the Hessian is a real symmetric matrix. Also, let the coordinate system for 
the space be defined by the orthonormal basis vectors {ei, . . . , e„}. 

We first choose the coordinate direction that is closest to the ridge normal 
direction. That is, we choose ej from the set {e^} such that |ej-vi| > 
|ej ■ vi| V 2 G {1, . . . ,n}. By Thm. Appendix A.l| we know that 

|ej ■ vi| > 1/y/n. 



(5) 



The gradient of F near xq may be written as 



1 



DF(x) = DF(xo) + ((x - xo) ■ D) DF(xo) + ^ ((x - xq) ■ D)^ DF(0 

for some ^ that is a linear combination of x and xq. Taking the dot product 
of DF(x) with ej and considering only points on a line normal to the ridge 
such that X = Xq + (ivi gives 

e,-DF(x) = e,-DF(xo)+t/e,-{(vi ■ D)DF(xo)} + ^rf2e,-{(vi ■ B^BFiO} • 

Since vi is an eigenvalue of Z)^-F(xo), this can be rewritten as 

e, ■ DF(x) = e, ■ DF(xo) + dX^e, ■ vi + ^d^e, ■ (vi ■ {D3F(0}) (6) 

We can establish an upper bound on ■DF(xo) by noting that vi ■DF(xo) = 
and expressing e^- and DF(xo) in terms of {vj}: 



|e,- ■ DF(xo) 



< 



< 



< 




5^(v, ■ DF(xo))v, 

i=2 
n 

5^(v, ■ DF(xo))v, 

i=2 

5^(v, ■ DF(xo))v, 



i=2 



^^(v. • DF(xo))^ 



i=2 



n 



|Di^(xo)||. 



(7) 
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The right hand side of Eq. [6] is a quadratic function in d with bounds on 
the coefficients 



< |c| < ci = VI - 1/n ||DF(xo) 



< 60 = ^ < |6| < 61 = |Ai| , 



n 



0<\a\< I -— < 

' ~ 4nv/^^^ DF xo 4ci 



Therefore, by Thm. Appendix A. 2 this function has a root at d* with 

^ 2v/^||DF(xo)|| ^ 

where ■ DF(x) = 0. 

We have now shown that the value of • DF(x) is zero somewhere on 
the fine normal to the ridge at xq at a distance of less than 

2v/(l-l/n) ||DF(xo)||/Ai. 

Call the zero point x*. To complete the proof, we will show that eJ(D^F)ej 
must be negative at this point so it is a coordinate local maximum. The first 
order Taylor series for the Hessian near xq is 

D2F(x) = D2f(xo) + ((x - xo) ■ D)D2F(0. 

At X* the second derivative of F in the ej direction is given by 

eJ(D2F(x*))e, = ej {B'F{^o))e, + rf*ej(vi ■ B'F{0)e, 

where d* < 1\/n — 1 ||DF(xo) || / Ai is the distance from the ridge of the point 
where ■ DF(x) = 0. Then 

e?D^F(x, < ^ + f 1 - i) A„ + ^^|lf<-°'" ||v. . D'F(0|| 

n \ nj |Ai| " " 

< 

where ||vi ■ D^F(^)|| is understood as the induced norm of the matrix that 
results from the tensor dot product vi ■ Y)^F{S,)- Thus, at x* the first deriva- 
tive of F is zero in the direction and the second derivative is negative, x* 
is a coordinate local maximum. 

□ 
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As long the criteria of Thm. 3.1 are satisfied along a ridge, there is a 
nearby coordinate local maximum that will be detected by the ridge tracking 
algorithm. The error of the location of the detected ridges is bounded by 



d 



2v/^^||DF(xo)|| 
|Ai| 



Since ||DF|| is typically 0{1) (if the ridge does not rise and fall too quickly), 
the error is typically C(l/ |Ai|). 

3.1. Ridge examples 

We now verify the properties of the ridge tracking algorithm on two two- 
dimensional examples. An analytically defined surface that admits a ridge 
and the FTLE field for a time dependent double gyre flow. We begin by 
examining the surface 



F{x,y) 



-0.1(x-j/)2 



ln((x + ?/)2 + 15)' 



(8) 



A contour plot of this function is 
shown in Fig. |4] and the function ad- 
mits a ridge along the line y = x. 
This represents a worst case scenario 
in terms of ridge orientation since the 
ridge is at a 45° angle to both coor- 
dinate directions. The ridge is also 
much thicker than is typically seen in 
FTLE ridges. Since this ridge is de- 
fined by an analytical function, it is 
possible to easily compute all the nec- 
essary criteria for Thm. 3^ In Fig. |4| 
the actual ridge has been drawn as a 
solid black line. The error bound, d is 
drawn in the figure as a pair of dashed 
black lines and the the coordinate local 
maxima have been drawn as red (for 
dF/dx = 0) and blue (for dF/dy = 0) 
lines. The other criteria of Thm 




Figure 4: The surface defined by Eq. [8] 
The solid black line denotes the ridge 
while the dashed black curves show the 
error bound and the red and blue curves 
show the coordinate local maxima. 



larger than is usually seen in FTLE ridges. We find Ai G 



0]are all satisfied despite Ai being much 

0.1477,-0.0843], 
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(a) 



(b) 



Figure 5: (a) The forward FTLE field and (b) the main FTLE ridge (with FTLE value 
> 0.25) for the time dependent double gyre with A = 0.1, e = 0.25, ui = 10, t — 0, and 
T = 15. 



A2 e [-0.0363,0.0071], and ||DF|| < 0.0356 on the ridge in the domain 

{x,y) E [-5,5]2. 

As a second example to apply the above theorem, we consider a ridge 
in the FTLE field of a time dependent double gyre. The velocity field con- 
sists of two counter rotating gyres with a periodic perturbation that enables 
transport between the two gyres. The flow is given by the stream function 



i'ix, y,t) = A sin(7r/(x, t))) sin(7r?/) 
on the domain [0,2] x [0, 1] where 

/(x,t) = a{t)x^ + b{t), 



(9) 



The velocity is given by 



a{t) 

m 



u = — 



e sin(ci;t), 

1 — 2e sin(a;t). 



(10) 



dy'' 



dil) 
dx 



(11) 



We use parameters A = 0.1, e = 0.25, and = 10 for this example and 
set the integration time for computing the FTLE at T = 15. We will only 
consider the forward time FTLE field at time t = which is shown in Fig. 

El 

The main FTLE ridge in this double gyre fiow (shown in Fig. |5}d) was 
computed with very high precision by iteratively estimating the ridge position 
and tangent direction and then adjusting the position in the normal direction. 
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Distance along the ridge Distance along the ridge 

(a) (b) 

Figure 6: (a) The FTLE values and gradient along the ridge shown in Fig. ^ and (b) 
bound on the ridge position error as given in Thm. |3.2| The distance along the ridge 
(x-axis) is computed as the distance along the ridge from the ridge origin at (1.110,0). 
Since ||VF|| = 0{1) and |Ai| ^ 1, the error bound is very small, typically less than 10~^. 



This is necessary to accurately compute the FTLE values on the ridge as well 
as the gradient, and Hessian of the FTLE field to bound the error in the ridge 
locations. Note that the ridge seen in this FTLE field is much sharper than 
the analytical ridge investigated in the previous example. 

Fig. [6] shows the FTLE values and the norm of the gradient of the FTLE 
field along the ridge as well as the bound on the ridge location error which 
is less than 10~^ for the majority of the ridge and never rises above 10~^. 
The norm of the gradient is bounded by ||V-F|| < L5. Additionally the 
eigenvalues of the Hessian fall in the range — L8 x 10^^ < Ai < —367 and 
-6.1 < As < 37.3 with averages of A^ = -8.1 x 10^° and A^ = -0.033 
and medians med(Ai) = —5.71 x 10^ and med(A2) = —0.30. The very large 
magnitude of Ai means that the grid based ridge tracking algorithm will be 
extremely accurate for this example. For reference, grid spacings of 10~^ or 
10~^ are typically used when performing LCS computations for this problem. 

4. Algorithm performance 

In this section we discuss the performance of the ridge tracking algorithm 
by examining two analytically defined examples: a time dependent double 
gyre and Arnold-Beltrami-Childress fiow. We verify that the expected LCS 
surfaces are extracted and then focus on establishing the computational order 
of the surface tracking algorithm and compare this to the standard FTLE 
algorithm that computes the FTLE field everywhere in the domain. 
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Forward (blue) and backward (red) LCS 




Figure 7: Surface tracking results for the double gyre flow. The 3D LCS are shown at in (a) 
with forward LCS colored blue and backward LCS colored red. The forward FTLE field 
is shown in (b) with the ridge tracking results overlaid as the black curve that precisely 
lines up with the FTLE ridge. 



4.1. Time dependent double gyre 

The first example presented is the time dependent double gyre used above 



in section We extend this fiow to three dimensions by simply setting the 
velocity in the z direction to w = 0. Since there is no z dependence and no 
velocity in the z direction, the LCS will be independent of z as well. We use 
the parameters A = 0.1, e = 0.1, and u = 27r/10. We set the integration 
time to be T = ±15 and compute both the forward and backward LCS. A 
threshold of 80% of the maximum FTLE value is used to determine the LCS. 

The LCS in this system have been well studied in the past and our results 
agree with previous publications [21 [H]- The full 3D LCS surfaces are shown 
in Fig. [7j This figure also shows the backward FTLE field overlaid with the 
results of the 3D ridge tracking algorithm. The LCS extracted by the ridge 
tracking algorithm lie exactly on top of the ridge in the FTLE field. 



The computational timing results are summarized below in Section 4.4 



4-2. Arnold- Beltrami- Childress flow 

Arnold-Beltrami-Childress (ABC) fiow is a three-dimensional, 27r-periodic 
fiow that has been previously studied with LCS techniques [13j. The fiow is 
given by 

u = A sm{z) + C cos(y), 

V = Bsm{x) + Acos{z), (12) 
w = C sm{y) + B cos(a;). 
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Forward (blue) and Backward FTLE field 




Figure 8: Ridge tracking results for ABC flow. The 3D LCS are shown in (a) with forward 
LCS colored blue and backward LCS colored red. (b) shows a single plane at z = tt to 
display the FTLE field (colored) and corresponding LCS (black curves) as computed with 
the ridge tracking algorithm. 

where A = I, B = .J^jl, and C = We use an integration time 

of T = 10 for FTLE computations and us a threshold value of 70% of the 
maximum FTLE value. 

The LCS are shown in Fig. |8] and agree with previously published re- 
sults [121 E] . Fig. |8] also shows single slice of the FTLE field well as the 
FTLE field and corresponding LCS at a height of z = it. The complex LCS 
present in the ABC flow are a product of the non-trivial invariant manifolds 
of this flow. They clearly divide the flow into different regions that appear as 
tube-like structures through the flow domain. These tubes are dynamically 
distinct from one another and particles travel within and along the tubes 
without escaping to other regions of the space. 

4 ■ 3. Swimming jellyfish 

As a final example, we compute the LCS created by a jetting type jellyfish. 
The jellyfish body motion was extracted from a digital video of a swimming 
individual jellyfish (species Sarsia tuhulosa) and use as input to an arbitrary 
Lagrangian-Eulerian CFD code that computes the resulting flow field and 
jellyfish acceleration. Full details of this procedure can be found in Sahin 
and Mohseni jl5j and Sahin et al. [16]. 

The resulting axisymmetric velocity field is returned on a moving, non- 
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uniform quadrilateral mesh in (r, 2;)-coordinates. During LCS computations, 
the {x, y, z) coordinates are converted to (r, z) coordinates to compute the 
velocity and the the velocity is converted back to Cartesian coordinates for 
particle advections. The mesh type creates significant complications for ve- 
locity interpolation during particle advection since even locating the mesh 
element that contains a given point is non-trivial. To address this issue ef- 
ficiently, an alternating digital tree (ADT) is used to search the domain for 
the element that contains a given drifter particle |17] . The ADT recursively 
divides the space in half so that at each node in the tree, only one branch 
must be searched. Since the nodes of the mesh elements are listed in counter 
clockwise order and the elements are convex, a point p is inside (including 
the boundary) an element if and only if 

Z ■ [(Vi - p) X (V(i n,od 4)+l - Vi)] > V i G {1, 2, 3, 4} 

for vertices {vj}. 

Once the element containing a particle is found, the velocity must be 
interpolated onto the drifter particle. Since the elements are generally none 
rectangular, simple linear interpolation is not possible. Instead, perspective 
projection is used to map the quadrilateral element and the point of interest 
onto the unit square. Bilinear interpolation is then used to approximate the 
velocity of the particle. 

The jellyfish chosen for this investigation has a swimming period of 1 s 
between contractions and an integration time of 0.5 s was used to compute 
the LCS. Fig. [9] shows the LCS as computed with the ridge tracking algo- 
rithm as well as the backward FTLE field computed with the standard FTLE 
algorithm. The backward LCS (Fig. |9^) clearly show a strong vortex being 
ejected as the jellyfish's bell contraction comes to an end. Additionally, the 
forward LCS outline fluid ahead of the vortex that will soon be entrained 
into the vortex ring as well as a region of fluid near the jellyfish bell that 
will be drawn into the bell during the relaxation phase of jellyfish swimming. 
These results are in excellent agreement with previously published LCS for 
this jellyfish |T8]. Furthermore, this example clearly demonstrates the visual- 
ization advantages offered by the ridge tracking algorithm. Since the actual 
LCS surface are computed it is much simpler to visualize both the forward 
and backward LCS simultaneously. 
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Forward (blue) and backward (red) LCS 




Backward FTLE field 




Figure 9: (a) The results of the ridge tracking algorithm and (b) the backward FTLE field 
for the swimming jellyfish. A strong vortex is being ejected near the end of the jellyfish's 
bell contraction. A cutaway view is shown so that the full LCS structure is visible. 
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4.4- Timing results 

We expect the computational time of the ridge tracking algorithm to be 
for a grid of spacing dx since computations are performed only 
near the 2D ridge surfaces. To establish the computational order, we have 
computed the full FTLE field using a standard FTLE algorithm and also 
computed the LCS with the ridge tracking algorithm presented above for the 
double gyre, ABC, and jellyfish flows. Computations were performed on at 
a variety of grid resolutions as well as on a single core and with a parallel 
code running on 16 or 48 cores. A least squares best fit was performed on a 
log-log scale for each case, assuming a fit of 

tf = C/{dx"). (13) 



The resulting data points and curve fits are shown in Fig. [10] and show that 
the standard algorithm is ~ 0{l/dx^'^) while the ridge tracking algorithm 
scales approximately as 0{l/dx'^-^). This performance is maintained for the 
parallel version of the code. 

It is worth noting that in the jellyfish example, for low resolutions the 
time required to read the velocity data files and build the ADTs used for 
search at each time step is a significant part of the total computational time. 
Since the velocity read in time does not change with grid resolution, this has 
the effect of giving artificially low exponents for the algorithm order (both for 
the standard algorithm and the ridge tracking algorithm). To compensate for 
this effect and more accurately estimate the asymptotic order, we estimated 
the velocity read in and ADT creation time and subtracted this value from 
the total run time before computing the computational order for the jellyfish 



example. These modified times are reported in Fig. 10 i. The resulting 
values for a closely match the values that would result from using only the 
last few data points and represent a closer approximation of the asymptotic 
values of a. 

It is also expected that the computational time should be directly pro- 
portional to the surface area of the LCS in the domain. We test this by 
generating an artificial FTLE field that has ridges along pre-selected planes. 
The planes are defined hj z = 0.1{x + y) + hi and the artificial FTLE field 
is given by 

a{x, y,z) = Y^ exp[-2000(2 - 0.1(a; + y) - hi)^]. (14) 
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ABC flow, 48 core ABC flow, 1 core 




10' 10'' 10' 10' 

1/dx 1/dx 

(a) (b) 



Double gyre, 48 core Jellyfish, 16 core 




10' 10' 10^^ 10^ 

1/dx 1/dx 



(c) (d) 

Figure 10: Timing results for the standard FTLE algorithm and the ridge tracking algo- 
rithm, a is the scaling exponent of the algorithm (CPU time « 0{l/dx°')) and appears as 
the slope of the lines in these log-log plots. The standard algorithm scales as 0{l/dx^ '^) 
and the ridge tracking algorithm scales as 0{l/dx^-^). 
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CPU time 


(s) 


Surface area 


1.01 


2.02 


3.03 


4.04 


5.05 


6.06 




1/64 


0.533 


0.689 


0.844 


0.987 


1.162 


1.333 


Grid Spacing 


1/128 


2.713 


3.306 


3.913 


4.625 


5.385 


6.121 


1/256 


19.54 


22.12 


25.10 


27.68 


32.79 


34.80 




1/512 


192.8 


204.5 


215.7 


235.7 


244.7 


259.9 



Table 1: CPU time (in seconds) to compute the artificial FTLE ridges of Eq. 14 for various 
surface areas and grid spacings. 



Each plane that defines an FTLE ridge is tilted slightly so that it is out of line 
with the grid to make detecting and tracking the ridge slightly more realistic. 
We use the domain [0, 1]'^ where each ridge has an area of Alcs = VU)2 for 
hi G (0, 0.8) and test six different cases corresponding to 1 — 6 ridges in 
the domain. Particles are advected using the double gyre velocity field listed 
above to account for the particle advection time, but instead of returning the 



true FTLE field for the double gyre fiow, the artificial field of Eq. 14 with 
pre-selected ridges is returned. The results are listed in Table [T] which shows 
the computational time according to surface area for four different values of 
dx. 

All four values of dx show a linear relationship between LCS surface 
area and CPU time and least squares fits result in the following regression 
coefficients for the fit tcpu = Ci + C2ALCS where Alcs is the surface area 
of the LCS and tcpu is the required CPU time: 



dx 


Ci 


C2 


1/64 


0.368 


0.157 


1/128 


1.94 


0.679 


1/256 


15.9 


3.14 


1/512 


177.9 


13.5 



The relatively large values of Ci for all these cases means that there is 
some initial cost regardless of the amount of LCS surface area. This is due to 
the cost associated with initializing data structures and performing the initial 



ridge detection step as described in Section |2.1[ Additionally, Ci appears to 
scale roughly as 0{l/dx^). Since the current implementation of the ridge 
tracking code uses and allocates full 3D arrays rather than using sparse data 
structures, it is reasonable to expect this relationship. On the other hand. 
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6*2 scales roughly as 0{l/dx'^). This accounts for the majority ridge tracking 
part of the algorithm. 

The 0{l/dx^'^) scaling of Ci may explain why the overall computational 
order of the algorithm is 0{1/ dx^'^) rather than 0{l/dx'^). As resolution is 
increased beyond current capabilities, the initialization cost of may be ex- 
pected to completely overwhelm the ridge tracking cost due to this difference 
in order. The implementation of sparse data structures would likely help 
solve this problem. 

5. Conclusions 

As the problems being analyzed with LCS techniques become increasingly 
complex, the corresponding computations become increasing expensive. The 
ridge tracking algorithm presented in this paper has been shown to reduce the 
order of LCS computations from 0{1/ dx^'^) to about 0{1/ dx"^'^) for three- 
dimensional flows. This reduction in order allows potentially tremendous 
savings in computational time as the required LCS resolution is increased. 

The effectiveness and algorithm properties have been demonstrated by 
several examples, including the analytically defined double gyre and ABC 
flow and the swimming jellyfish that is defined by velocity stored in data 
files. On single processor as well as multicore machines, the ridge tracking 
algorithm shows the expected change in computational order and provides 
large speed ups for all tested cases. 

We have also proved that although the ridge tracking algorithm detects 
coordinate local maxima rather than the actual ridges, for well defined ridges 
the associated error is small. The distance between a ridge and the surfaces 
detected by this algorithm is (9(||DF|| / |Ai|) where DF is the gradient along 
the ridge and | Ai | is the second derivative normal to the ridge (or, equivalently 
the smallest eigenvalue of the Hessian of F). In typical examples this error is 
at smaller than the grid spacing and for well defined ridges it may be several 
orders of magnitude smaller than the grid spacing. 

The general framework of this algorithm could easily be adapted to other 
LCS definitions or techniques. For example, the finite size Lyapunov expo- 
nent (FSLE) could easily be used in lieu of the FTLE. A 2D version of this 
algorithm (or the grid-less algorithm of Lipinski and Mohseni [TT]) could be 
relatively easily modified to use the variational formulation for LCS which 
was recently proposed by Haller [6l [19] . 
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Future work will focus on implementing sparse data structures to further 

reduce the computational cost of large LCS calculations and on implement- 
ing some addition LCS techniques such as those mentioned in the previous 
paragraph. We expect that the use of sparse data structures may further 
reduce the order of computation from 0{l/dx'^-^) nearer to 
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Appendix A. Additional proofs 

This appendix contains additional theorems and proofs referenced by the 
proofs of the theorems in the text. 
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Theorem Appendix A.l. Given an orthonormal basis {ej} for M" and 

an arbitrary unit vector v there is a basis vector ej that maximizes |v • ej| 
and |v • ej| > Ij ^/n. 

Proof. Existence of the maximum is implied by the extreme value theorem 
since the dot product is a continuous function and the set {sj} is finite and 
therefore closed and bounded. The other part of the proof is as follows: 
Assume |v • e^l < Then 





n 
i=l 






n n 

(E(ei • v)ei) • ■ v)ei) 

1=1 i=l 




n 

i=i 






n 

\ i=i 




< V 


/n(l/x/^)2 




< 1 







which is a contradiction so |v • ej| > 1/y/n. □ 

Theorem Appendix A. 2. Given a quadratic equation y — ax^ -\-bx -\- c 
with bounds on the real valued constants 

|c| < Ci 
< 6o < l&l < bi 

4ci 



\a\< '° 



there is a real root x* such that \x*\ < 



6, 







27 



Proof. If a = 0, there is a single root at x* — —c/b and = \—c/b\ < 
2ci/6o. 

If a 7^ 0, the roots are given by the alternate form of the quadratic formula 



2c 

X — 



b±Vb^-iac 
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The discriminant, A = 6^ — 4ac > 6q — 4— ^ci = is positive so the quadratic 

4ci 

has two real roots. Denote these roots 



Xl = 



X2 



\x — 



2c 


-b + 


- Aac 


2c 




-b-Vy^- 


- Aac 


t X* — X2- 


Then 


2|c| 




\b\-V^ 


- 4ac 



2ci 
^~b^ 



□ 
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